Derivatives of the fixed-node energy 



O 

o 

(N 

> 
O 



> 
oo 



O 



o 



S. De Palo('^'^), S. Moroni^"^'') and S. Barom(°''=) 
"■Istituto Nazionale per la Fisica della Materia 
Dipartimento di Fisica, 
Universita di Roma "La Sapienza" 
'^Scuola Internazionale Superiore di Studi Avanzati, Trieste 

Italy 

We present quantum Monte Carlo calculations of total energy derivatives, consistently performed 
in the fixed-node approximation. Contributions from nodal displacements, neglected or approxi- 
mated in previous investigations, are properly taken into account. Their impact on the efficiency is 
discussed. 



I. INTRODUCTION 

Quantum Monte Carlo simulations are emerging as a 
promising alternative route towards high-accuracy pre- 
dictions of molecular_and materials properties. The fixed- 
node approximationld is widely used in quantum Monte 
Carlo simulaljonsa of many-fermion systems to avoid the 
sign problenjj. It consists of solving for the lowest-energy 
eigenstate, $, of the Hamiltonian, H, with the boundary 
conditions that the wave function vanishes on the nodal 
surface of an anti-symmetric trial function, ^, i.e. at the 
points R in configuration space where '^{R) = 0. The 
fixed- node energy (FNE) — i.e. the eigenvalue E of the 
Hamiltonian corresponding to $ — is an upper bound for 
the ground state energy, and it becomes exact if \1/ has 
the same nodes as the true ground state wave function, 
^ex- Therefore, for a given trial function, the FNE is 
more accurate than the variational energy (the expecta- 
tion value of H on ^P), because exactness of the latter 
requires the stronger condition that coincide with ^ex 
for all R. Indeed, the fixed-node approximation offers the 
best viabie numerical solution for several many-fermion 
problemgj. 

The ability to calculate low-order derivatives of the 
Born-Oppenheimer, energy with respect to an external 
parameter. A, is an important ingredient of the calcula- 
tion of materials properties. Forces acting on individ- 
ual atoms or the stress acting within an extended sys- 
tem are examples of first-order derivatives whose vanish- 
ing determines the mechanical equilibrium of a system. 
Similarly, vibrational frequencies — as well as any kind 
of static generalized susceptibility — are determined by 
second-order derivatives. The calculation of derivatives 
of the FNE, dE/dX, has the potential of largely extend- 
ing the scope of applications of quantum Monte Carlo 
simulations and spawning further methodological devel- 
opments (besides the examples given above, A could also 
be the strength of any external field, or a variational pa- 
rameter which modi£es the nodal surface of ^'). However 
technical difficulticso have long hindered routine calc 
tion of FNE derivatives, and despite recent progress 
open questions still remain. In particular, the nodal dis- 
placement upourvariation of A has always been neglected 
in toto or in partB, thus adding a further bias on top of the 



fixed-node approximation: referring to recent literature, 
A-independent nodes have been chosen in the path inte- 
gral Monte Carlo calculation of electronic forceaS; they 
have to be choserujao in applications of the Hellmann- 
Feynman theoremQEI; finally, the drift-diffusion factor of 
the propagator, together with its nodal contributions, has 
been dropped from the weights m correlated-sampling 
quantum Monte Carlo simulationaj. 

In this paper we address the calculation of the deriva- 
tives of the FNE without any further approximations. 
To this end, we examine in detail a simple but significant 
model that, we believe, is representative of the difficulties 
which are met in the calculation of these derivatives in 
more realistic situations. The focus being on the role of 
nodal displacements, we consider a system where such a 
role is crucial and a check against exact results is possi- 
ble, namely a spinless Fermi gas in two dimensions in a 
weak external field. 



II. METHODOLOGICAL ASPECTS AND THE 
MODEL 



The Hamiltonian of our model system is: 
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= -;^-Vf + A^cos(qT,), 
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where the strength of the external field. A, is the pa- 
rameter upon which the FNE is considered to depend, 
and where periodic boundary conditions have been as- 
sumed. The ground state is a Slater determinant, D, 
of A-dependent Mathieu functionstj, whereas our trial 
function is = J x D with a symmetric Jastrow fac- 
tor J(i?) = exp(— J u{rij)) which introduces spurious 
pair correlations. The fixed-node solution, exact in this 
case, is recovered as an absolutely non-trivial result of the 
simulation, since the trial function, wrong on purpose, 
can be made poor at will by tuning the pseudo-potential, 
u. 

We use| the reptation quantum Monte Carlo (RQMC) 
method", in which the calculation of derivatives is sim- 
pler to implement than with branching algorithms — just 
as simple as in path integral Monte CarlcO. The FNE 
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can be computed from the mixed estimate of the Hamil- 
tonian, 



E 



mH\<^) _ jEL{R)f{R)dR 



jJ{R)dR 



(2) 



by averaging the local energy El{R) — H'^{R)/'ii{R) 
over the mixed distribution /(i?) = <i>(i?)^'(i?). The un- 
known state $ is expressed in terms of 5* through the 



imaginary-time propagator, $ = lim,-- 



^^f. Upon 



discretization of the imaginary time evolution and intro- 
duction of importance sampling, the mixed distribution 
can be cast into the form: 



/(i?M)= lim /n^oiG(i?„i?,+i;e)*(i?o)' 

M-^OO J 



(3) 



where e — t/M is the time step. The 
importance-sampled Green's function, G{R, R';e) = 
^'(i?)(i?|e-"^|i?')/*(i?'), is replaced Jpy a short-time 
approximation, a common choice beingB 

GfIR, R'- e) oc e-(«'-«-^^(«)/2)V2.g-f 

(4) 

with F = Vlnf^, 

In the RQMC method a generalized Metropolis al- 
gorithm is used to sample the probability distribution 
for a path of finite length in imaginary time (a reptile, 
X — {i?o, ■ • ■ , Rm}, 

P{X) = n^-iG(i?„ it'.+i; e)*(i?o)', (5) 

the hmit in Eq. (||) being reached for large enough M . 
Replacement of Eq. (||) in (||) shows that the FNE can 
be evaluated by averaging the local energy on the distri- 
bution P, 

E {{euRm))) = jdxp{x) • 

Since P{X) is explicitly known, we can differentiate to 

get 

dxE = {{A{X))), (7) 
A{X) = [dxEL{RM) + {EL{RM)-E)dxliiP{X)] 

Finite increments of amplitude A are used for the evalu- 
ation of d\EL{R) and 9a \nP{X) on the path configura- 
tions sampled from P{X). Equivalently, finite increments 
can be used for E via the correlated-sampling method: 



^, ^ JdXP'{X)Ei{RM) 
JdXP'iX) 

^ J dXP{X)WiX)E'^iRM) 
JdXP{X)W{X) 



(8) 



with W{X) = P'{X)/P{X), primed quantities being 
evaluated at A + A. 

In the variational Quantum Monte Carlo (VMC) ap- 
proximation, the probability distribution P{X) in Eq. (0) 




FIG. 1: Importance-sampled Green's function for a one- 
dimensional motion near a hard wall. Left panel: exact; mid- 
dle panel: Eq. right panel: Eq. with a smoothed F. 
The range of d and d' is (0, e^/^). 



or (||) is replaced by ^'^{Rm)- Small energy differences 
are routinely calculated by VMCu. The more accurate 
FNA requires instead the full probability distribution of 
Eq (|^), which entails nontrivial difficulties through the 
term dx In P{X) of Eq. (|^) . Therefore we now examine in 
some detail the behavior of the probability distribution P 
when a configuration along the path, say Ri , approaches 
a node of 

When the distance d from the nodal surface is small 
enough, the dynamics of the system is dominated by 
the hard \ra,ll potential due to the fixed-node boundary 
conditionaij. The importance-sampled Green's function 
should behave like that of a one-dimensional motion close 



to a hard wall [Go(c^, d'; e) cx -je 



§!_p-(d-d' f 111 



(1- 



as shown in the left panel of Fig. (Q) . The short time ap- 
proximation Gf{R, R' instead, is a very poor approx- 
imation near a node, because it assumes F = V In to 
be constant during a time step, whereas F diverges as 
d~^. The approximation (Q) for a one-dimensional mo- 
tion close to a hard wall, shown in the middle panel of 
Fig. (|l]), also bears little resemblance to Gq. First, Gp 
does not vanish at ci = where Go does. Sampling from 
Gp would thus imply an unphysical finite density at the 
nodes; the correct density, which is proportional to d^ 
can be recovered including in the algorithm a Metropolis 
test based on the square of the trial functiontfEJ. Sec- 
ond, Gf vanishes at d' = where Go does not. This 
would imply a divergence proportional to d^^ of InP at 
the nodes, thus making the evaluation of the derivative, 
Eq. (|^), impossible. Such a divergence can be eliminated, 
leaving the zero-time-step limit unchanged, by a suitable 
cutoff on the quantum force, Vln^. We thus replace 
F in Eq. | with F = 4F(-1 + yTTF^e/2) / {F^ ejB. 
The resulting Green's function, Gp , is shown in the right 
panel of Fig. (nh. With this replacement the variance of 
the estimator (H) of the derivative is finite. However we 
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empirically find that the extrapolation of the result at 
e — > can be problematic (see also Ref. [10]) 

In the simulation of electrons in atoms and molecules, 
it proves convenient to use an interpolation between 
Eq. (01 and a locally better Green's function near the 
nucleili3. In the same spirit, we use an interpolation G 
between Gp and a locally better Green's function Gn 
near the nodes: 

G{R, R';e) = a{d)Gp{R,R';e) + (l - a{d))GN{R,R';e) 

(9) 

where the function a{d) smoothly varies from to 1 as c? 
goes from to ~ e^^^. In the new Green's function 

Gj,{R,R';e) oc |gle(«'-«)V4e(i _ ,-.'.A) (iq) 

the potential V is treated in the primitive approximation, 
whereas the nodal dependence has the correct behavior 
of Gq. Although the variance due to d\ \nP{X) is finite 
using G, a time-step dependent cutoff on d is legitimate 
and can be useful. 



TABLE I: First derivative of fixed-node energy calculated 
with Variational Monte Carlo (VMC), using RQMC with G 
and using the approximated algorithm as explained in the 
text 

Exact VMC RQMC App. Alg. wave-function 
1 0.910(16) 0.998(38) 0.941(14) Wi 
1 0.9677(92) 1.038(38) 0.997(13) *// 
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III. RESULTS 

We turn now to the discussion of our results for the 
Fermi gas. We first choose a very poor trial function ^P/ 
with a pseudo-potential, u = aexp(— 6r^) where a = 0.9 
and b = 0.7, yielding a variational energy about 20% 
higher than the exact ground state energy (whereas the 
FNE and the exact result coincide in this case, since 
the Slater determinant is the true ground state wave 
function). The derivative dxE\x(,, with Aq = 0.02 and 
q — 1.585, is calculated using G. The result, shown in 
Tab. ^ does indeed recover the correct value. 

For comparison, we ako perform calculations using the 
approximation of Ref. which drops the drift-diffusion 
part of the propagator in Eq. (||) and introduces an ex- 
tra factor ^'^(i?j\f) in order to recover the exact result 
in the limit of perfect importance sampling. As shown 
in Tab. |, the approximate algorithm is more efficient (it 
gives a statistical error 3 times smaller for the same num- 
ber of path configurations). However its result for the 
derivative is clearly biased. Using a better trial function 
^'// with a = 0.45 the systematic bias of the approx- 
imate algorithm becomes smaller than statistical error. 
It is interesting to note that the ratio of the statistical 
error between the two methods does not seem to depend 
strongly on the quality of 5*. We also have evidence that 
the efficiency depends on the number of particles only 
through the extra cost of sampling configurations for a 
larger system (simulations have been performed for 5, 9 
and 21 particles). 

In Fig. ^ we show the time-step dependence of vari- 
ous derivatives calculated using different Green's func- 
tions. Data obtained with Gp have smaller statistical 



FIG. 2; Time step dependence of the first derivative of the 
FNE {dxE\x=o.oo2), upper panel, and of the density fluctu- 
ation (9ApqU=o)) lower panel. We show calculations using 
G (black squares and diamonds) and Gp (gray circles and 
stars); the black lines are quadratic fits, whereas the gray 
lines include a square root term. 

errors than with G, but their extrapolation to zero time- 
step is harder to guess in advance. The upper panel 
shows derivatives of the FNE, 9a£^|a=o.oo2- The calcu- 
lation with G extrapolates to the exact result assuming 
a quadratic time-step dependence, whereas an e^^^ term 
seems necessary in order to account for the Gp data. 

The problems with the e — > limit using Gp are 
more evident for the derivative of the density fluctua- 
tion, dxpqlx^o, where Pq(-R) = I]iexp(-iq- r^): as we 
can see from the lower panel of Fig. (g), in this case an 
even smaller power of e would be required to extrapolate 
to the exact result at e = 0. 

The physical interest in d\Pq\x=Q stems from its rela- 
tion to the static dielectric response. The derivative of 
Pq is evaluated as {{B{X))), with_B(A) = 9pq(i?A//2) + 
Pq{RM/2)d\\nP{X) (cfr. Eq. ^jB). The same physical 
information can be formally obtained from d^E\\^Q. Sec- 
ond derivatives are in principle straightforward to evalu- 
ate, but differentiation of Eq. (^) leads to terms contain- 
ing 9^ InP and {d\ InP)^ which have large fluctuations. 

Despite its low efficiency, the calculation of second 
derivatives of the energy offers an interesting conceptual 
advantage over the first derivative of the density fluctua- 
tion. Suppose we don't know the exact nodal dependence 
of 'i' on the strength A of the external field. Starting from 
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FIG. 3: Second derivative of the FNE, dxE\{^2o ^ function 
of a. The parameter of the pseudo-potential is a = 0.90, and 
the calculation uses the Green's function G. The exact result 
corresponds to a = 1, d^E]^^^^^ = —1. 



plane- wave nodes at A = 0, we can explore a family of 
wave functions, in which the one-particle orbitals are 
Mathieu functions calculated for an external potential of 



strength aX different from the physical value A. The 
nodal dependence of a on A varies with a: for a = 0, 
no nodal displacements are allowed; for a — 1, the exact 
nodal dependence is recovered. We can compute within 
the same run (hence in a correlated way) the values of 

^"x^txlo pertaining to several choices of a. The result is 
shown in Fig. (^). Because of the variational character 
of the fixed node approximation, the optimal nodal de- 
pendence (within the family ^'q,, and within error bars) 
is given by the value a corresponding to the minimum 
of d\E\^^2f^. Inspection of Fig. (^ leads us to conclude 
that: (i) suppression of the nodal dependence (a = 0) 
leads to a small statistical error, but produces a large 
bias; (ii) the minimum of S^-EI^'^'q is consistent with the 
exact nodal dependence {a — 1); (iii) the dependence of 

^a-^Ia=o quadratic, so that calculations for just 

3 values of a would suffice (this could be useful for a 
higher-dimensional optimization) . 

Unfortunately, the calculation of second derivatives 
has still very large statistical errors (if nodal displace- 
ment is properly included) . We hope that ajSuccessful im- 
plementation of noise reduction techniquestS could even- 
tually turn it into a convenient computational tool. Fur- 
ther perspectives in nodal optimization are offered by the 
calculation of first derivatives of the FNE, for instance in 
conjunqtipn with the stochastic gradient approximation 
methodllZl. 
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